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Metastatic carcinoma cells exhibit at least two different phenotypes of motUity and invasion - amoeboid and 
mesenchymal. This plasticity poses a major clinical challenge for treating metastasis, while its underlying 
mechanisms remain enigmatic. Transitions between these phenotypes are mediated by the Racl/RhoA 
circuit that responds to external signals such as HGF/SF via c-MET pathway. Using detailed modeling of 
GTPase-based regulation to study the Racl/RhoA circuit's dynamics, we found that it can operate as a 
three-way switch. We propose to associate the circuit's three possible states to the amoeboid, mesenchymal 
and amoeboid/mesenchymal hybrid phenotype. In particular, we investigated the range of existence of, and 
the transition between, the three states (phenotypes) in response to Grb2 and Gabl - two downstream 
adaptors of c-MET. The results help to explain the regulation of metastatic cells by c-MET pathway and 
hence can contribute to the assessment of possible clinical interventions. 



Metastasis is a hallmark of cancer', and responsible for more than 90% of solid tumor deaths". Carcinoma 
is the most common type of solid tumor, whose metastasis is a complex process that begins when some 
epithelial cells from the primary tumor undergo epithelial-to-mesenchymal transition (EMT) to lose 
cell-cell adhesion and gain migratory and invasive mesenchymal characteristics. Carcinoma cells can deploy 
different strategies (phenotypes) to tread through the extra- cellular matrix (ECM)'', and different interactions 
with the stroma or local microenvironment' (Fig. la). Our previous theoretical investigations revealed that the 
core regulatory circuit of EMT operates as a three-way switch, allowing not only for epithelial (E) and mesench- 
ymal (M) phenotypes but also for a hybrid epithelial/mesenchymal phenotype (E/M), which is associated with 
collective cell migration*'\ Yet, the plasticity of the individual cell selection between different motility character- 
istics (phenotypes) stLU remained illusive. It is widely accepted that deciphering the underlying mechanisms of 
cellular plasticity during metastatic invasion is central for designing therapeutic targeting of carcinomas'. To help 
meet this challenge, we present here theoretical investigations of the GTPase-based operation principles of the 
Racl/RhoA circuit - the key regulator for amoeboid-to-mesenchymal transition (AMT). The modeling challenge 
is to simplify the complexity of this circuit whose dynamics involves transcription, translation, and post- 
translation (GTPase) regulations. 

Phenotypic Plasticity of Individual Cell Migration. Carcinoma cells typically adopt two different phenotypes to 
invade the three-dimensional (3D) matrix environment - the mesenchymal (M) or the amoeboid (A) (here we 
specifically refer to blebby amoeboid (BA))'' (Fig. la). Cells of the M phenotype are elongated and spindle-shaped 
with their leading edge characterized by lamellopodia (LAM) and/or filopodia (FIL). These cells are able to 
remodel and even degrade the ECM by secreting Matrix Metalloproteinases (MMPs) and therefore act as 'path 
generators''''^. Conversely, rather than secreting MMPs to remodel the ECM, cells that exhibit the A phenotype 
have higher shape plasticity that enables them to squeeze into the gaps in the ECM, thus acting as 'path finders''"''. 
In 3D environment, many carcinoma cells exhibit mesenchymal-to-amoeboid transition (MAT) and amoeboid- 
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Figure 1 | Plasticity of cell migration phenotypes and the core regulatory circuit, (a) The phase diagram showing the Relationship between the activities 
oftheGTPases and the plasticity of ceU migration phenotypes. Background colors correspond to the different level of Racl-GTP and RhoA-GTP - red for 
high level of Racl-GTP, green for high level of RhoA-GTP, yellow for high levels of both of them, and blue for low levels of both of them. The activity of the 
GTPases is hypothetically associated with the cell morphology and mobility. The phenotypes are depicted as cartoons displaying their corresponding 
morphological features - epithelial phenotype (E), hybrid epithelial-mesenchymal phenotype (E/M), mesenchymal phenotype (M), amoeboid phenotype 
(A), and hybrid amoeboid-mesenchymal phenotype (A/M). M phenotype is characterized with lamellopodia (LAM) and/or filopodia (FIL), and A 
phenotype here is specifically referred to Blebby amoeboid (BA) phenotype, which is characterized by blebbing. The A/M phenotype is considered as a set 
of different morphologies - Lamellipoida with blebs (LB), Lobopodia (LP) and Pseudopodal amoeboid (PA). The blue color is also associated with strong 
cell-cell adhesion, as observed in E or E/M phenotypes, while the rest colors are associated with single cell migration modes, (b) The core AMT/MAT 
regulatory circuit connected with the c-MET pathway. The AMT/MAT is mainly regulated by the Rac 1 /RhoA regulatory circuit, while RhoA and Rac 1 are 
regulated via Grb2 and Gabl by the c-MET pathway, which receives the external signal from HGF/SF. A solid arrow denotes activation, and a solid bar 
stands for repression. A solid line represents transcriptional regulation, and a dashed line is for non-transcriptional regulations such as GTP loading, (c) 
Dynamical system characteristics of the Racl /RhoA regulatory circuit. The plot shows the nuUclines and possible steady states corresponding to Equation 
(1). Without any external signal (Grb2 = 0, Gabl = 0), the circuit can be tristable. Rednullclineisford[i{c*]/A = Oandblacknullclineisford[i{(,']/A = 0. 
Green solid circles denote the stable fixed points, and green hoUow circles denote the unstable fixed points. Each stable point can be associated with a cell 
phenotype, depicted as a cartoon beside them. 
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to-mesenchymal transition (AMT) either spontaneously*'* or in 
response to external signals from the microenvironment*' '"'". 

In recent studies, carcinoma cells have been observed to have 
phenotypes with hybrid amoeboid/mesenchymal (A/M) character- 
istics (Fig. la)*'"* '^"'. For instance, some carcinosarcoma cells have 
both lamellopodia/filopodia and blebs structures (LB)**'"", which have 
also been seen in vivo during early development of the developing 
Fundus" and zebrafish embryos"*. Some fibroblasts adopt cylin- 
drical-shaped lobopodia phenotype (LP) in the 3D environment 
with morphological features of both lateral blebs and blunt protru- 
sions", and some leukocytes and neutrophils show pseudopodal 
amoeboid (PA) migration phenotype in 3D as well as 2D environ- 
ments with dynamic protrusions in the front part and high contract- 
ility in the rear part" '\ Such a rich phenotypic plasticity during the 
migration of individual cell enables the tumor cells to adapt to their 
changing microenvironments, and plays a crucial role in cancer 
dissemination". 

The Racl/RhoA GTPase-based Regulatory Circuit. The choice 
among the aforementioned phenotypes is operated by the Racl/ 
RhoA regulatory circuit. Racl and RhoA belong to the Rho family 
of small GTPases and act as molecular switches by changing between 
their active (the GTP-bound state) and inactive (the GDP-bound 
state and the GDI-bound state) forms^". This switching process is 
regulated by three sets of proteins: GEFs (Guanine Nucleotide 
Exchange Factors) that catalyze the exchange of bound GDP for 
GTP, thus elevating the levels of the active GTPases; GAPs 
(GTPase-activating proteins) that promote the intrinsic GTP 
hydrolysis rates, thus reducing the concentration of the active 
GTPases; and GDIs (Guanine nucleotide dissociation inhibitors) 
that sequester GTPases from the membrane to the cytosol and 
stabilize the proteins by preventing degradation^'. 

Racl and RhoA regulate the phenotypic transitions by controlling 
the actin polymerization and actomyosin contraction"'^. Therefore, 
the activities of these two GTPases have been observed to correlate 
with cell morphology and motility. For example, actomyosin con- 
tractility increases in response to the RhoA activation, thus resulting 
in membrane blebbing and facilitating the amoeboid phenotype 
^y^j9,23.24 jj^g other hand, the activation of Racl results in the 
formation of focal adhesions and actin polymerization, which leads 
to the formation of lamellopodia and enables a mesenchymal pheno- 
type (M)"'"'' "'^. Appropriate changes in the relative strengths of these 
two driving forces (actomyosin contraction vs. actin polymerization) 
allow for not only the transitions between the A and M phenotypes, 
but also might enable the transition into/from the hybrid A/M 
phenotype**. These transitions can be triggered by extracellular sig- 
nals such as Hepatocyte Growth Factor/Scatter Factor (HGF/SF) 
(Fig. lb) through the c-MET pathway^'. c-MET, the specific tyrosine 
kinase receptor for HGF/SF, is often overexpressed in many carci- 
nomas and correlates with poor patient survivaP''^". The activated c- 
MET recruits Grb2 (Growth-factor receptor bound protein 2) and 
Gabl (Grb2 associated binding protein 1) to regulate the activity 
of both RhoA and Racl (See Supplementary Table SI and 
Supplementary Fig. SI for the detailed molecular interactions). It is 
contradictory that the HGF/SF/c-MET pathway has been reported to 
be able to induce either mesenchymaP'"*" or amoeboid"" phenotype, 
therefore the underlying mechanisms of the HGF/SF/c-MET path- 
way remain elusive. 

Importantly, the active GTP-bound forms of Racl (RhoA) inhibit 
the activation of RhoA (Racl) by promoting the hydrolysis of the 
active GTP-bound forms of RhoA (Racl) (See Supplementary Table 
SI), thus acting as a mutually inhibitory feedback circuit, often 
named as a toggle switch. The Racl/RhoA feedback circuit may have 
unique dynamical properties, because the regulation of the circuit is 
post-translational in nature, which is distinct from the dynamics of 



either transcriptional (TF-TF) or transcriptional and translational 
(miR-TF) regulation in well-studied canonical toggle switches*'^ '^''^. 

While the RhoA-Racl circuit has been extensively studied experi- 
mentally' "'"'' it has received limited theoretical attention, prim- 
arily due to the special small GTPase-based regulation. Here, we 
developed a new framework to model the small GTPase-based regu- 
lation in the context of the Racl/RhoA circuit by considering detailed 
transitions among different states of GTPases. We note that this 
article considers the expression and overall activity level of the 
GTPases. More specifically, in this first step model, we did not 
include spatial effects of the GTPases - the fact that Racl and 
RhoA are activated in different spatial compartments of the cell. 
Inclusion of the spatial effects is important and will be done in future 
extension of the current model. The model simulations revealed that 
the Racl -RhoA circuit can act, for a wide range of realistic para- 
meters, as a three-way (ternary) switch between three possible states: 
the amoeboid phenotype (A), the mesenchymal phenotype (M) and 
the hybrid amoeboid/mesenchymal phenotype (A/M). The model 
demonstrated that Grb2 and Gabl could induce the activation of 
both Racl and RhoA, which is expected to promote the migration 
of cells through mesenchymal and amoeboid phenotypes respect- 
ively. Since both Grb2 and Gabl are regulated by the c-MET path- 
way, this may possibly explain the observation in which HGF/SF/ 
c-MET pathway induces either of these two mutually exclusive phe- 
notypes in different experiments"^'". We will discuss several experi- 
mental observations that are consistent with our new model. The 
model provides the first step towards understanding how different 
levels of the small GTPases, RhoA and Racl, in a cell govern pheno- 
typic plasticity during carcinoma metastasis under the influence of 
external signals in the tumor microenvironment. 

General View of the Core Regulatory Unit. For clarity, we wiU first 
summarize the key findings of our model on Racl/RhoA circuit, 
followed by more detailed presentation of the analysis. In Fig. la, 
we have shown a spectrum of motility/ invasion phenotypes during 
carcinoma cell metastasis. The epithelial (E) and epithelial/ 
mesenchymal hybrid (E/M) phenotypes are not discussed in this 
article, since we focus on the motility stage that the cells have 
already undergone complete EMT. The core AMT/MAT regulatory 
unit, as shown in Fig. lb, is composed of mutual inhibitions between 
the active forms of RhoA and Racl (RhoA-GTP and Racl-GTP). This 
mutually inhibitory circuit is driven by signals from the c-MET 
pathway through Gabl and Grb2 (See Supplementary Table SI). 

We first considered the stand-alone dynamics of the Racl/RhoA 
circuit. We showed that the circuit acts as a three-way switch for a 
wide range of biologically realistic parameters (Fig. Ic) similar to the 
characteristic behavior of other 'self- activating toggle switches' 
(SATS)"*^. In the absence of any external signal, the three states of 
this circuit (Fig. Ic) are: 1. (1, 0) state - high RhoA-GTP with low 
Racl-GTP that corresponds to A phenotype, as characterized by 
blebs (BA). Notably, in some cases, high RhoA-GTP is reported to 
inhibit cell migration depending on ECM stiffness and cell type''^'^''; 2. 
(0 + , 1) state - low RhoA-GTP with high Racl-GTP that corresponds 
to M phenotype as identified by the presence of lamellopodia (LAM) 
and/or filopodia (FIL). Here, we denoted the level of RhoA-GTP as 
"0 + " to take into account the fact that mesenchymal cells still need 
some active RhoA in the rear part for retraction", while active Racl is 
present at the leading edge. Therefore, in mesenchymal phenotype, 
the expression levels of the active forms of RhoA and Racl are spa- 
tially separated'", whereas in amoeboid phenotype, active RhoA is 
almost uniformly distributed''^; 3. (1, 1) state - balanced relatively 
high RhoA-GTP and Racl-GTP that we propose to be consistent 
with the hybrid A/M phenotypes observed experimentally"''''''^ "'"""'. 
The existence of the hybrid (1, 1) state indicates richer migration 
plasticity. The predicted multistability in this circuit, i.e. the existence 
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Figure 2 | Schematic diagram for theoretical framework and its appUcation on Racl-RhoA regulatory circuit, (a) Schematic diagram of the regulation 
of a typical Rho family GTPase (denoted as Rho). The inactive GDP-bound state of Rho (Rho-GDP) and the active GTP-bound state of Rho (Rho-GTP) 
both bind to the membrane. They can interconvert through the regulations of GAPs (at rate /[/J) or GEFs (at rate B[/2])> which maybe activated by some 
external input signals (Ii and I2). Rho-GDP can be released from the membrane by binding to a GDI molecule (at rate g(ii_R) and revert to its membrane- 
bound state by releasing GDI (at rate dgdi_R). Rho-GDP and Rho-GTP degrade at rate Kj^ and Kr,, respectively, while the degradation of Rho-GDI was not 
considered, because GDI binding can stabilize the Rho protein^', (b) The RhoA-Racl regulatory circuit (Details in Supplementary Table SI). The GTP- 
bound states of RhoA or Racl can promote GTP loading of its own, and meanwhile activate the GTP hydrolysis of the other. RhoA-GTP is also 
transcriptionally self-activated. Grb2 induces the GTP loading of Racl, while Gabl induces that of both Racl and RhoA. (c) The effective (reduced Racl/ 
RhoA) circuit. In terms of Racl-GTP and Rho-GTP, their mutual inhibitions form a non-canonical toggle switch with positive auto-regulations. A solid 
arrow denotes activation, a solid bar indicates repression, and the wavy line represents regulation by external signals. The solid double line represents the 
GTP loading or hydrolysis process while the dashed double one represents the binding or unbinding process of GDI molecules. The dashed-dot lines 
indicate the indirect regulations on GTP loading or hydrolysis process via GEFs or GAPs. (d) Typical values of the B and the j functions with respect to the 
concentrations of the GTPases. The B and J functions represent the GTP loading and hydrolysis rates of both Racl and RhoA respectively. Both functions 
increase with the level of GTP-bound Racl and RhoA. The parameters for the two functions were listed in Supplementary Table S2. 



of diverse possible states, can explain the observed existence of the 
diverse phenotypes during cell migration. 

Effective Model of the Small GTPase-based Regulatory Circuit. 

The challenge posed in modeling the small GTPase-based 
regulatory circuit is to incorporate the elaborate transitions 
between different forms of the GTPases. Typically, a GTPase 
protein can switch among its active (GTP-bound state) and 
inactive (GDP-bound state and GDI-bound state) forms under the 
regulation of three sets of proteins (GEFs, GAPs and GDIs)''". To 
understand these features of the GTPases, we developed the 
theoretical framework for GTPase-based circuit by specifically 
modeling the detailed molecular interactions as illustrated in 
Fig. 2a. We utilized the approach to investigate the dynamics of 
the core Racl/RhoA regulatory circuit, as shown in Fig. 2b. 

The deterministic dynamics of the circuit could be modeled by a 
set of six chemical rate equations presented in the Supplementary 
Equations (7) and (8). Yet under the assumption that the total level of 
Racl or RhoA (the sum of levels of GTP-bound, GDP-bound and 
GDI-bounded form) always reaches a steady state, the above detailed 
model can be approximated by an effective model (Fig. 2c) described 



by the following two rate equations (See more detail in Supplementary 
Information Section 2): 

at Kd^ 

{H[K])+f{[K])+h„+Ku^)[K], 

where Rc* represents the active Racl (Racl-GTP), and represents 
the active RhoA (RhoA-GTP). is the production rate for Racl, gi^i^ 
and gKi^A are basal and excitatory production rates for RhoA respect- 
ively. Kr> and Krj are their corresponding degradation rates. Ir^ and 
Iri^ represent two external signals that drive the circuit, such as Grb2 
and Gabl in c-MET signaling. The HiU function^ ([_R/,*]) represents 
the transcriptional self- activation of RhoA. Two new functions, B and 
J, are defined to represent the GTP loading and hydrolysis rates 
(See details in Supplementary Information Section 2). The values of 
the both fimctions monotonicaUy increase with an increase in the 
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concentration of the active GTPases, as illustrated in Fig. 2d. The 
effective model was used for stability and bifurcation analysis while 
the detailed model was used for the dynamic simulations. The model 
derivation and the parameter estimations are described in details in 
the Supplementary Information Section 3. 

The Racl/RhoA Circuit as a Three-way Switch. We started by 
analyzing the circuit dynamics in the absence of external signals 
(Grb2 =0 and Gabl = 0). As illustrated by a typical phase-plane 
diagram in Fig. Ic, Racl/RhoA regulatory circuit can act as a three- 
way switch among the following three states: (high active RhoA/low 
active Racl), (low active RhoA/high active Racl), and (both high 
active RhoA and Racl), which we denote as (1, 0), (0 + , 1) and (1, 
1) respectively. According to the experimental observations"'^^"^^'''", 
we associate the states (1, 0) and (0 + , 1) with the amoeboid (A) and 
mesenchymal (M) phenotypes. Again, we used "0 + " to denote some 
minimal level of active RhoA present in the rear end of mesenchymal 
cells and required for their individual migration^^. Furthermore, we 
proposed to associate the (1, 1) state with the amoeboid- 
mesenchymal (A/M) hybrid phenotype that has been suggested in 
some recent experiments both for cancer cells and normal cells 
during early embryonic development"*''''^'"''^""*. Although these 
experiments lack quantitative measurement of the activity of Racl 
and RhoA, the mixed morphologies of these cells share the traits of 
both amoeboid and mesenchymal phenotypes. These properties may 
be indicative of relative high levels of both active RhoA and Racl . Yet, 
it is clear that a direct measurement of these proteins is necessary to 
establish a direct association between these phenotypes and the 
model predicted hybrid phenotype. 

Inhibition 
External Signal 



Notably, for a small range of parameters, the Racl/RhoA circuit 
can also act as a four- way switch, which gives rise to the coexistence 
of four states or quadra- stability (Supplementary Fig. S2a). The 
fourth state could be associated with the E/M hybrid phenotype that 
exhibits collective cell migration because activation of Racl and 
RhoA are needed for EMT^''''^ and balanced, intermediate levels of 
Racl and RhoA activity are suggested in experiments on the E/M 
hybrid phenotype or partial EMT'*'' Since we mainly focused on 
solitary movement, we limited our analysis on the parameter range in 
which Racl/RhoA circuit acts as a three-way switch. 

The Switch Response to External Activation and Inhibition 
Signals. Next, we analyzed the response of the Racl/RhoA circuit 
to an external input signal that drives the circuit through either Racl 
or RhoA. We modeled that the signal directly increases or decreases 
the basal OTP loading rates. In Fig. 3, we show the response to the 
external signal that affects the RhoA loading rate by the effective 
model shown in Fig. 2c. When the signal activates the basal GTP 
loading rate for RhoA, it gives rise to the coexistence of the diverse 
phenotypes A, M, A/M, all of which correspond to solitary 
movement. However, when the signal inhibits the loading rate, it 
also gives rise to collective cell migration of the E/M hybrid 
phenotype. As is also evident from Fig. 3, a high activation of the 
GTP loading rate leads the cell to a monostable phase {A} in which 
only amoeboid phenotype (A) exists, whereas low levels of GTP 
loading rate correspond to the monostable phase {M} in which 
only the mesenchymal phenotype (M) exists. The activity level of 
Racl for each phenotype is shown in Supplementary Fig. S3a. Our 
model is consistent with the experiments showing that cells with 

Activation 
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0.1 0.2 0.3 0.4 0.5 
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Figure 3 | Bifurcation of RhoA-GTP protein levels in response to an external signal regulating the GTP loading rate of RhoA. The signal can 
either increase (activation) or decrease (inhibition) the GTP loading rate. The response to activation/inhibition is shown on the right/left side of the 
bifurcation respectively. The red solid lines indicate stable states and the blue dashed lines indicate unstable states. The bifurcation illustrates the possible 
coexistence (for some range of the signal) of four states: (i) the (1, 0) state with high RhoA-GTP and low Racl-GTP, which corresponds to A phenotype; 
(ii) the (0 + , 1), which corresponds to M phenotype; (Hi) the (1, 1), which correspond to A/M phenotype; (iV) the (0, 0) state, which corresponds to E/M 
phenotype. The corresponding bifurcation of Racl-GTP protein levels is shown in Supplementary Information Section 5. Co-existence of different 
phenotypes is associated with a multistable phase, highlighted by different background colors (see legend at the bottom). Starting from the (0-I-, 1) state 
(M phenotype, at bottom left part of the red curve), the system stays in the M phenotype when the inhibition signal is reduced; further switching the 
inhibition signal to an increasing activation signal leads the system to undergo a transition to the (1, 0) state (A phenotype, at top right part of the red 
curve). The transition is indicated by the dashed upward arrow at the boundary of the phase {A, M} and {A}. Similarly, increasing the inhibition signal can 
induce the transition from the (1,0) state (A phenotype) back to the (0 + , 1) state (M phenotype), as indicated by the downward arrow at the boundary of 
the phase {M} and {A, M}. Besides, cells may switch to the A/M or E/M phenotype due to noise in gene expression. 
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constitutively active RhoA are associated with amoeboid (A) or 
blebby amoeboid (BA) phenotype'"; while cells with dominant 
negative RhoA usually exhibit a mesenchymal (M) phenotype'". It 
also predicted that, when the external signal acts only on RhoA, the 
induced transitions between A and M phenotypes are much easier 
than the transitions from A or M phenotype to A/M or E/M 
phenotype. This may explain why these hybrid phenotypes, A/M 
and E/M, were rarely observed during AMT in most of the 
experiments". 

We have shown the response to an external signal connected to 
Racl (see Supplementary Information Fig. S3b,c). In this case, inhibi- 
tion of the Racl loading rate did not give rise to the existence of the E/ 
M hybrid phenotype. However, this phenotype does exist when two 
external signals of inhibiting nature drive both Racl and RhoA 
simultaneously (Supplementary Fig. S3d). These results indicate that 
the cells may attain the E/M phenotype when GTP loading rates of 
RhoA or Racl are low (see blue area in Fig. la). However, when these 
GTP loading rates increase significantly due to the external signals, 
the cells are more likely to switch to one of the solitary modes of 
migration thus displaying reduced plasticity of spontaneous transi- 
tions among different migration phenotypes, which is consistent 
with experiments from Turner group"*^ (see green, red and yellow 
areas in Fig. la). 

The Switch Response to Input Signals from Grb2 and Gabl. As 

was mentioned earlier, Grb2 activates only Racl, and Gabl activates 
both Racl and RhoA. To understand the circuit response to these 
regulations, we first investigated the response of the circuit dynamics 
to either Grb2 or Gabl (in terms of the corresponding bifurcation 
diagram) when they act individually. We found that when Grb2 level 
is increased, the cells adopt a mesenchymal (M) phenotype 
(Supplementary Fig. S4a,b); whereas Gabl induces the cell to 



adopt an amoeboid (A) phenotype (Supplementary Fig. S4c,d). 
However, further high Gabl signal can induce the cell to adopt the 
amoeboid/mesenchymal (A/M) phenotype since Gabl can activate 
both Racl and RhoA (Supplementary Fig. S4e,f). 

Next, to understand the combined effect of Grb2 and Gabl, we 
constructed the two-parameter bifurcation diagram, as shown in 
Fig. 4a. Each phase corresponds to a particular situation in which 
one or several different phenotypes can coexist. More specifically, the 
possible phases are: 1. Phases with only one phenotype - {A}, {M} 
and {A/M}. 2. Phases in which two phenotypes can coexist - {A, A/ 
M}, {M, A/M} and {A, M}. 3. A phase in which all three phenotypes 
coexistence - {A, M, A/M}. The detailed information of the various 
phases indicates the plasticity of cell migration as driven by different 
combinations of Grb2 and Gabl signals. Depending on how Grb2 
and Gabl increase temporally, the cells follow different trajectories in 
this phase diagram and thus go through different phenotypic transi- 
tions as is illustrated in Fig. 4b and 4c. 

The regulatory effect of Gabl depends on the relative strengths of 
its activation to the GTP loading rates of Racl and RhoA (details are 
presented in the Supplementary Section 9). In the case of stronger 
activation of RhoA, Gabl stimulates the cells to acquire the amoeb- 
oid phenotype. In contrast, in the case of stronger activation of Racl, 
Gabl stimulates the cells towards the A/M hybrid phenotype or the 
mesenchymal phenotype (Supplementary Fig. S7). Understanding 
the circuit response to input signal from Grb2 and Gabl provides 
crucial clues regarding the pleiotropic effects of the c-MET pathway 
in promoting either the amoeboid or mesenchymal mode of migra- 
tion and also transitions between them (AMT and MAT). 

Phenotypic Transitions Driven by the c-MET pathway. To better 
understand the phenotypic transitions, we investigated the response 
dynamics of the Racl /RhoA circuit when the input signals, Gabl and 
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Figure 4 | The circuit response to input signals from Grb2 and Gabl. (a) Two-parameter bifurcation phase diagram. Two input signals Grb2 
(x-axis) and Gabl (y-axis) are selected as the two control parameters. As explained in the text, eacli phase corresponds to a different combination of 
coexisting plienotypes (Phase plane diagrams for each phase are shown in Supplementary Fig. S5). For example, the orange area is the phase {M, A/M}, 
which means that cells in this phase can belong to either the M or the A/M phenotype. The colors used for different phases are explained by the legend in 
Fig. 3. (b) One-parameter bifurcation diagram for the circuit driven by Gabl when the Grb2 level is fixed to be 8 X 10"^ molecules. The corresponding 
trajectory in the two-parameter bifurcation phase diagram (Fig. 4a) is shown by a purple line. The transitions between the M and A/M phenotypes are 
illustrated by the dashed upward and downward arrows, (c) One-parameter bifurcation diagram for the circuit driven by Grb2 when the Gabl level is 
frxed to be 2.5 X 10"^ molecules. The corresponding trajectory in the two-parameter bifurcation phase diagram (Fig. 4a) is shown as a brown line. The 
transitions between the A and A/M phenotypes are illustrated by the upward and downward arrows. 



SCIENTIFIC REPORTS | 4 : 6449 | DOI: 1 0. 1 038/srep06449 



6 





4 200 

I 
I 
I 



400 



0 0.2 0.4 0.6 0.8 1.0 
Grb2(x 10' molecules) 

•'fx 




S 2.1 



0 0.2 0.4 0.6 0.8 1.0 
Grb2(>!l0' molecules) 



O 



600 800 idoo 1200 
Time (Hours) 




1 






lA 


1 A/M j 


M 




1^ N 




-Racl-GTP [ 1^ 






1 


-RhoA-GTP 1 1 






/ 









4.0 
.1.2 
2.4 
1,6 
0.8 



0 200 400 600 800 1000 1200 
Time (Hours) 



Figure 5 | Temporal dynamics of the cells to HGF/SF treatment, (a) The time dependent Gabl and the Grb2 signals that imitate the effect of 
HGF/SF treatment, as explained in the text, (b) Trajectory of the cell response projected on the phase diagram. The results are for a cell in which Gabl 
activation of RhoA GTP loading is stronger than its activation of the Racl GTP loading (The parameters are gtp_Ri,l2 = 240/i"' and gtp_Rcl2 = 90h^' 
respectively). The solid line is the trajectory, and both the arrows on the line and the color gradient of the line (from black to white) indicate the time 
evolution. The blue star marks the initial condition, the green star marks the peak of Gab 1 expression and the pink star marks the peak of Grb2 expression, 
(c) Time dynamics of the expression levels of Racl-GTP (Red) and RhoA-GTP (Black) in response to the signals, Gabl and Grb2. The different 
phenotypes during the transition are highlighted by different background colors, where green, yellow and red areas stand for A, A/M and Ivl phenotype 
respectively, (d) Similar to (b) but for a cell in which Gabl activation of Racl loading is stronger than its activation of the RhoA loading (The parameters 
are gtpJR^h = lOOO/i"' and gtp_Ri,l2 = 240^"' respectively), (e) Similar to (c) but for the case shown in (d). 



Grb2, change in time. To recapitulate the possible response of the 
circuit to HGF/SF treatment, we chose time dependent functions for 
the levels of the Gabl and Grb2 to mimic the cells' response to HGF/ 
SF treatment (Fig. 5a). The HGF/SF signal leads to increase c-MET 
phosphorylation"^'', which in turn regulates Gabl and Grb2 together 
with Met-Induced Mitochondrial Protein (Mimp) in a form of two 
coupled feed-forward loops (FFLs) as is shown in Fig. Ib'"* '''^ Hence, 
the levels of Gabl and Grb2 were modeled as two pulse signals with a 
time delay as shown in Fig. Sa""". 

The form of the pulses is similar to that of a typical response of a 
FFL. Time delay is incorporated to reflect the effect of the feed- 
forward like coupling of Grb2 and Gabl to c-MET (Gabl responses 
ahead of Grb2). The simulated treatment starts with cell in the 
amoeboid (A) phenotype at the left bottom corner (shown as blue 
star) in the phase diagram shown in Fig. 5b. The cell stays in this 
phenotype when Gabl is increased (Fig. 5c) but makes a transition 
into the hybrid (A/M) phenotype after Gabl decreases and Grb2 
increases, thus the activity of Racl increases whUe that of RhoA 
remains almost unchanged. Finally, after Grb2 also decreases, the 
cell goes through another transition from the hybrid (A/M) pheno- 
type into the mesenchymal (M) one. The results illustrate how c- 
MET pathway can regulate the cells to switch between different 
migrating phenotypes. In Fig. 5d and 5e, we demonstrated the result 



for the same simulation but for a ceU with different circuit parameters 
(a cell in which Gabl activation of Racl GTP loading is stronger than 
its activation of the RhoA loading). In this case, the increase of Gabl 
signal induces the cell to transit into the hybrid A/M phenotype 
instead of maintaining it in A phenotype. The different signal- 
response behaviors of these two simulations support that different 
cell lines (cells with different parameters) may respond differently to 
the regulatory signals Gabl and Grb2. 

The aforementioned results are consistent with experimental 
results and help to explain several experimental observations. In 
particular, c-MET pathway is reported to induce both the mesench- 
ymal phenotype for several non-cancer and cancer cells""'''™ and the 
amoeboid phenotype for some breast cancer cells^^'. Other experi- 
ments show that Grb2 is essential for TGF-P to induce a mesench- 
ymal phenotype for some cancer cells^' and Gabl can stimulate AMT 
by forming dorsal ruffles through its adaptor Nck^^. These observa- 
tions have been well captured by our simulations showing that 
mesenchymal phenotype can be induced either by Grb2 in one cell 
line or by Gabl in another cell line with different parameters (see 
Fig. 5 and Supplementary Fig. S7). 

Testable Predictions. We present here our predictions that could be 
tested in future experiments. Input from Grb2 stimulates the cells 
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towards mesenchymal phenotype (Fig. 6a) while input from Gabl 
stimulates the cells towards the amoeboid phenotype (Fig. 6b). These 
predictions can be tested, in principle, by changing the expression of 
Grb2 and Gabl in given cells. For a certain set of parameters, which 
can be considered representing some cell lines in which Gabl has 
stronger activation of Racl GTP loading, Gabl is predicted to 
stimulate the cells to become mesenchymal or the hybrid phenotype. 

Phenotype Distribution. Cells belonging to the same cell line often 
display non-heritable phenotypic variabUity''^^'''^. Such variability can 
originate, for example, from local differences in the microenviron- 
ment (such as ECM rigidity) leading to differences in the circuit 
parameters of the individual cells. We have shown earlier that cells 
with different circuit parameters can respond in a different way to the 
input signals. Hence, we expect to see a distribution of phenotype for 
given input signals. As a first step to assess the expected nature of the 
population level distribution, we extended our simulations to a 
population of 5,000 cells, each with different circuit parameters. 
More specifically, the cells that compose the populations have 
±5% variations from the original parameters. 

In Fig. 7, we showed the percentages of cells that can be in one of 
the three different possible phenotypes, (A, A/M and M), for differ- 
ent levels of the input signals. We found that for high Grb2 or high 
Gabl signal, a significant percentage of cells adopt the A/M pheno- 
type (Fig. 7a,b). This result, which is obtained due to a weak robust- 
ness (high sensitivity to the circuit parameters), indicates that cells 
under high Grb2 or Gabl signal are still sensitive to the conditions of 
their microenvironment. However, both high Grb2 and Gabl signals 
can result in more cells being in the hybrid A/M phenotype (Fig. 7c). 
We also generated cell phase distribution for different initial condi- 
tions and found similar results (See Supplementary Fig. S8). 

Discussion 

By considering the relationship between active levels of Rho GTPases 
and cell morphology, we were able to build a theoretical framework 
on Racl/RhoA GTPase-based regulatory circuit to interpret some 
experimental observations about cancer cell migration phenotypes, 
and further present several testable predictions for future experiments. 



In the current model, we devised the circuit to capture the GTP 
loading and hydrolysis based mechanism of mutual inhibition 
between the GTPases Racl and Rho A. As is showed in Supple- 
mentary Table SI, most of the connections have been discovered 
in breast cancer cells, but a few of them are found in the other cell 
lines. It is important to note that in some cell types, the circuit might 
be different; for instance, some positive feedback between Racl and 
RhoA is also reported in 3T3 fibroblasts via Dbs and mDia. However, 
this positive feedback seems to be cell-type specific and the fun- 
damental mechanism stiU remains unknown"'". Therefore, we 
assumed that the circuit we devised can be relevant to many types 
of cancers, and thereby can provide valuable understanding of the 
underlying mechanism of AMT in general. 

We found the Racl/RhoA regulatory circuit is usually a three-way 
switch. Such tristability has been shown to be a hallmark of many self- 
activating toggle switches (SATS), i.e. double negative feedback loops 
with self- activation on one or both of the elements of the loop''"''''^^'^''. 
Here, based on experimental observations''^'"^^''"', we propose to 
associate the (1, 0) - high activity of RhoA but low activity of Racl 
state with the amoeboid phenotype, and to associate the (0+, 1) state 
- low activity of RhoA and high activity of Racl with the mesench- 
ymal phenotype. Notably, it has been reported that high activity of 
RhoA can also inhibit cell migration by forming strong focal adhe- 
sions to ECM, that is specific to some microenvironments'^ and cell 
types'"". In the current model, we did not consider the interplay 
between cells and ECM. Hence, it is not possible to distinguish the 
migratory and stationary phenotypes in the case of high RhoA activ- 
ity within the current theoretical framework presented here. The 
interaction between the Racl -RhoA circuit and the ECM as mediated 
via integrins remains a subject for the future extension of the model. 

In addition to the (1,0) state (A phenotype) and the (0 + , 1) state 
(M phenotype), the model predicts the existence of a third (1,1) state 
with balanced relative high RhoA-GTP and Racl -GTP. We proposed 
to associate this state with the experimentally observed"'^'^''''''^""' 
hybrid ameboid/mesenchymal (A/M) phenotype. As was explained 
earlier, the reason for this association is as follows: both high RhoA- 
GTP and Racl-GTP can enable cells in this state to obtain morpho- 
logical properties resulting from both high actomyosin contractility 
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Figure 7 | Phenotype distribution of a population of cells at different 
levels of Grb2 and Gabl signals, (a) Phenotype distribution of a population 
of cells driven by Grb2 signal (Gabl = 0). The cell parameters are randomly 
distributed over ±5% relative to the original parameters. The y-axis denotes 
the percentage of cells corresponding to a specific phenotype. The color of 
each line represents a phenotype. Under the stimulation of different level of 
Grb2, a population of cells has different phenotype distribution (percentage 
of cells shown in y-axis). For each point, the original phase of the cells 
depends on the level of Grb2 signal. Note that due to the co-existence of 
different phenotypes in one phase (e.g. in both the phases {A, M, A/M} and 
{A, M}; A and M phenotypes are present), the sum of the total percentages 
for one particular signal level can be more than 100%. For example, the 
initial phase with neither Gabl nor Grb2 signals is {A, A/M, M}. About 
100% cells can be A or M phenotype, and about 36% cells can be A/M. Some 
of these cells were in some multistable regime comprised of one or more of 
the phenotypes of A, M and A/M. (See Supplementary Fig. S8). (b) 
Phenotype distribution driven by Gabl signal (Grb2 = 0). (c) Phenotype 
distribution under Grb2 and Gabl regulations. The color of each bar, 
corresponding to the color definition above represents the initial cell phase. 
When both Grb2 and Gabl signals are high, the cells are highly likely to be 
maintained in hybrid A/M phenotype. 

and high actin polymerization, such as pseudopodal amoeboid 
phenotype mentioned above. Moreover, the downregulation of the 
effect of Racl or RhoA could induce the transition from A/M to 
canonical A or M phenotypes" '^"". For example, the knockdown of 
RhoA or inactivation of myosin II switched cell migration from 
lobopodia (LP) or filopodia/blebs (LB) to exclusive lameUipodia 
(M) mode"'". Since the (1, 1) state for different cell types might 
manifest itself as different phenotypes'* ''''''''' '''"', we lumped them 
together as A/M phenotypes in this work. 

Regarding implications for cancer cell migration, we note that the 
cells belonging to the hybrid phenotype are expected to be more 



adaptable when cells migrate in a complex environment. More spe- 
cifically, since these cells can display some properties of both amoeb- 
oid and mesenchymal phenotypes, and can also promote the 
transitions between these two phenotypes, they can more readily 
adopt to ECM with high variability. The model prediction presented 
here call for future direct measurements of both active levels of RhoA 
and Racl in these phenotypes to further prove our hypothetic asso- 
ciation of the (1,1) state with these A/M phenotypes. 

We also showed the diverse regulatory functions of Grb2 and 
Gabl on Racl /RhoA regulatory circuit. Since Grb2 and Gabl are 
two different adaptors for c-MET receptor, their different functions 
in inducing cell migration may help explain how c-MET pathway can 
induce the cells to adopt either of these A or M phenotype, which are 
usually mutually exclusive. It may be noted that Grb2 and Gabl both 
form Feed-Forward Loops (FFLs) with c-MET receptor, thus their 
activity may be separated in time. In Fig. 5, we show that AMT can be 
induced by sequential activation of Gabl and Grb2 signals to the 
regulatory circuit. This indicates that integrating c-MET pathway in 
detail with Racl /RhoA regulatory feedback loop in future might 
better explain the impact of HGF/SF on cell migration. 

In this study, we have explored in detail the phenotypic plasticity 
pertaining to single cell migration. However, when sessile epithelial 
cells from primary carcinoma undergo EMT, they start moving col- 
lectively in the hybrid E/M state, and some of them further transit to 
move individually in mesenchymal or amoeboid mode; therefore 
either completing EMT""'^ or undergoing a Collective to Amoeboid 
Transition (CAT)". The dynamic regulation between E, E/M, A, M 
and A/M phenotypes remains far from understood. Future investi- 
gations into understanding the coupled dynamics of EMT/MET and 
MAT/ AMT regulatory networks - miR-200/ZEB and RhoA-Racl 
respectively - would be imperative in charting out the entire plas- 
ticity landscape that migrating carcinoma cells can adopt"". 

As noted before, the current model does not incorporate intracel- 
lular diffusion of the GTPases, i.e. the effect of spatial distributions of 
these Rho GTPases. Since these effects can be significant, an import- 
ant future extension of the model should be to incorporate the intra- 
cellular spatial dynamics of Racl and RhoA and its connection with 
the cell shape change between the three different phenotypes - A, M 
and A/M. More specifically, this would entail understanding the 
spatial distribution of these proteins within cells exhibiting different 
phenotypes, for instance, cells in M phenotype has high levels (1) of 
active Racl on the cell side extending the protrusion, and some 
minimal amount of active RhoA (0 + ) on the opposite pole, and 
the diffusion of molecules between different ends of the cell may 
introduce a time delay. The role of Cdc42, a Rho GTPase that is 
important for Racl localization and the formation of filopodia""^', 
would also be decisive in this endeavor to understand how spatial 
distribution of activity levels of RhoA and Racl govern the cell mor- 
phology. Previous models for spatial segregation of Rho GTPases, 
based on reaction-diffusion mechanism, have explored the role of 
mutual inhibition and fast diffusion in stabilizing the cell shape^""™. A 
recent attempt in this direction used a simple Boolean model of 
RhoA and Racl activity in two different spatial compartments in 
the ceU'". The Boolean framework considers the activity levels of 
RhoA and Racl to be binary - either ON (1) or OFF (0), and does 
not capture the existence of any intermediate state(s). However, we 
show that the Racl /RhoA loop is a three-way switch that allows for 
the existence of hybrid A/M phenotypes. Therefore, our framework 
is better equipped to chart out comprehensively the shape space of a 
cell based on active levels of Racl and RhoA. 

To conclude, our model provides a better understanding of the 
plasticity of cell migration and its regulation by external signals. 
Furthermore, it paves a promising way to understand how c-MET 
pathway is involved in carcinoma metastasis. Given the large 
number of current attempts to therapeutically target the c-MET 
pathway, understanding this relationship can provide some useful 
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